2.1 Dataset for the volcano plot

Dataset we’ll use comes from an experiment where mice were exposed to x-ray (this is also used as a treatment to kill cancer cells) . Then messenger RNA (mRNA) levels were measured from ~27,000 genes in spleen tissue in both control and exposed mice. The table below shows the results from a statistical test comparing these two groups across ~27,000 genes. Conducting the same comparison test with 27,00 genes leads to (multiple comparison (testing) problem)[https://en.wikipedia.org/wiki/Multiple_comparisons_problem]. To overcome this problem and detect the true positives, multiple testing correction was applied to the results, which gives the false discovery rate (FDR).

genes.Spleen.Xray <- read_csv("data/results_miceSpleen_xray.csv")
head(genes.Spleen.Xray)

This data format is already in appropriate format for our plotting purposes.

We will only add columns:
- logFDR: Negative \(log_{10}\) transformed FDR (false discovery rate) values
- threshold: defining whether a gene qualifies for the set threshold. In this case, it is any gene that satisfies:
- \(FDR < 0.1\) (which is also same as \(logFDR > 1\)) and \(|\ logFC\ | > 0.4\)

We will also pick the top 10 genes and assign as a separate data frame for labels in the plot.

pdat <- genes.Spleen.Xray %>%
  mutate(logFDR = -log10(FDR),
         threshold = ifelse(FDR < 0.1 & abs(logFC) > .4, TRUE, FALSE) )

#subset top 10 genes to add labels later to our plot
data.label <- pdat %>%
  filter(threshold == TRUE) %>%
  top_n(n = 10, wt = logFDR)
data.label

2.2 geom_point layers

In the plot any point below our threshold has different aesthetics mapped to them!! Therefore, we will add two layers of geom_point and provide different datasets and aesthetic mappings to each layer.

One of the differences between these two layers is the shape of the point.

Plotting ‘character’ (pch) There are multiple shapes available for plotting points. Take a look with:

?pch

Back to our plot…

vp <- ggplot(pdat, aes(x = logFC, y = logFDR)) +
  geom_point(data = filter(pdat, threshold == FALSE),
             size=.35, alpha = .55, colour = "grey80") +
  geom_point(data = filter(pdat, threshold == TRUE),
             aes(size = logFDR, fill=logFDR),
             pch = 21, stroke = .25, alpha = .65, color = "grey59")
vp

2.3 Let’s talk about colors

?brewer_pal
brewer.pal(6, "PuRd")

## [1] "#F1EEF6" "#D4B9DA" "#C994C7" "#DF65B0" "#DD1C77" "#980043"

display.brewer.pal(6, "PuRd")

myPallete <- colorRampPalette(brewer.pal(8,"PuRd"))
myPallete(10)

##  [1] "#F7F4F9" "#EAE5F1" "#DCCAE3" "#D0ACD3" "#CB8EC4" "#DC6AB2" "#E43C96"
##  [8] "#DB1E72" "#C00E50" "#91003F"

show_col(myPallete(10))

show_col(myPallete(100), labels = FALSE)

2.4 Customize color scales, legend & axes labels and plot title

vp <- vp +
  scale_fill_gradientn(guide="colourbar", colours = myPallete(100)) +
  guides(size=FALSE,
         color=FALSE,
         fill=guide_colourbar(title = expression(-log[10](FDR)))) +
  labs(x = "Fold Change",
       y = expression(-log[10](FDR))) +
  geom_hline(yintercept=1, color="black", size=1) +
  annotate(geom = "text", x=-6, y=0.75, label="italic(FDR) == 0.1",parse=TRUE) +
    ggtitle("A Volcano Plot!")

vp

2.5 Add labels for top 10 genes

vp <- vp  +
  geom_label_repel(data = data.label,
                   aes(label = GENE), alpha = .95,
                   size = 4, color = "black",
                   box.padding = 0.5, segment.size = .5,
                   arrow = arrow(length = unit(0.015, 'npc')))


vp

2.6 Make cosmetic changes with theme()

vp <- vp  +
  theme_minimal() +
  theme(plot.title = element_text(lineheight = 1,face = "bold",size = 20),
           axis.text = element_text(size = 12),
           axis.title = element_text(size = 16, face = "bold"),
           legend.text=element_text(size=13),
           legend.position="bottom",
           legend.title = element_text(size = 12,face = "bold"),
           legend.key.size = unit(0.035, "npc"),
           strip.text.x = element_text(size = 16, colour = "black"),
           strip.text.y = element_text(size = 16, colour = "black"))

vp

2.7 Save plot

Default is PDF format but many graphic devices are available.

?ggsave

ggsave( "volcano_plot.pdf",
        vp, width = 10, height = 8)

ggsave( "volcano_plot.jpeg",
        vp, width = 10, height = 8, device = "jpeg")

3.1 mpg data

ggplot2 dataset, you need to first load tidyverse or ggplot2

Fuel economy data from 1999 and 2008 for 38 popular models of car

?mpg
head(mpg,3)

mpg$class[1:20]

##  [1] "compact" "compact" "compact" "compact" "compact" "compact" "compact"
##  [8] "compact" "compact" "compact" "compact" "compact" "compact" "compact"
## [15] "compact" "midsize" "midsize" "midsize" "suv"     "suv"

3.2 Barplots

geom_bar is designed to make it easy to create bar charts that show counts. For example number of cars in each class

ggplot(mpg, aes(manufacturer)) +
  geom_bar() +
    theme(axis.text.x = element_text(size = 10, angle = 45))

ggplot(mpg, aes(manufacturer)) +
 geom_bar(aes(fill = class)) +
  theme(axis.text.x = element_text(size = 10, angle = 45))

ggplot(mpg, aes(manufacturer)) +
 geom_bar(aes(fill = class), position = "dodge")

ggplot(mpg, aes(manufacturer)) +
  geom_bar(aes(fill = class),position = position_dodge(width = 0.4))

3.3 Visualizing data points or mean as barplot

Plotting mean or other statistics of a variable directly with ggplot2 is less straightforward. One way is to use stat_summary_bin function.

?stat_summary_bin

ggplot(mpg, aes(class)) +
  stat_summary(aes(y = displ), fun.y = "mean", geom = "bar")

A more straightforward approach is to use geom_col.

df <- data.frame(trt = c("a", "b", "c"), outcome = c(2.3, 1.9, 3.2))
df

ggplot(df, aes(trt, outcome)) +
  geom_col()

ggplot(df, aes(trt, outcome)) +
  geom_col(fill="steelblue")

ggplot(df, aes(trt, outcome)) +
  geom_col(aes(fill=trt))

mpg %>%
  group_by(class) %>%
  summarise(cty_mean = mean(cty)) %>%
  ungroup() %>%
  ggplot(aes(class, cty_mean)) +
  geom_col()

3.4 geom_bar() with continuous data

You can use geom_bar() with continuous data, in which case it will values and show counts at unique locations. (But geom_histogram is more appropriate for this task.)

ggplot(mpg, aes(displ)) + geom_bar()

5.1 Example ToothGrowth dataset

head(ToothGrowth)

?ToothGrowth

Here’s how the the distribution of tooth length for different regimens look like.

gb <- ToothGrowth %>%
      unite(regimen, -len, sep= "_") %>%
      ggplot(data=., aes(regimen, len)) +
      geom_boxplot(aes(color=regimen))

6.1 Compute stats for ToothGrowth

tg.stats <- ToothGrowth %>%
  group_by(supp, dose) %>% # group by
  summarise(.,
            mean.len = mean(len),  # mean
            sd.len = sd(len),      # sd
            n = length(len)        # count
            ) %>% 
  ungroup() %>%
  mutate(se.len = sd.len/sqrt(n))

Now this is how tg.stats look

tg.stats

Finally we can start plotting! Let’s start with line graphs.

6.2 Point and Line graphs

Step by step:

ggplot(tg.stats, aes(x=dose, y=mean.len, colour=supp)) +
    geom_point(aes(x=dose, y= mean.len), size=4, shape=18)

ggplot(tg.stats, aes(x=dose, y=mean.len, colour=supp)) +
  geom_point(aes(x=dose, y= mean.len), size=4, shape=18) +
  geom_errorbar(aes(ymin=mean.len-sd.len, ymax=mean.len+sd.len), width=.1)

ggplot(tg.stats, aes(x=dose, y=mean.len, colour=supp)) +
  geom_point(aes(x=dose, y= mean.len), size=4, shape=18) +
  geom_errorbar(aes(ymin=mean.len-sd.len, ymax=mean.len+sd.len), width=.1) +
  geom_line()

ggplot(tg.stats, aes(x=dose, y=mean.len, colour=supp)) +
  geom_point(aes(x=dose, y= mean.len), size=4, shape=18) +
  geom_errorbar(aes(ymin=mean.len-sd.len, ymax=mean.len+sd.len), width=.1) +
  geom_line() +
  # here I add the original data points!
  geom_point(data=ToothGrowth, aes(x=dose, y=len, fill=supp), alpha=0.4, size=1.5)

7.1 1. Change titles and set scale limits

Modify axis, legend, and plot labels see Reference ggplot2: Scale

ggplot(tg.stats, aes(x=dose, y=mean.len, colour=supp)) +
  geom_point(aes(x=dose, y= mean.len), size=4, shape=18) +
  geom_errorbar(aes(ymin=mean.len-sd.len, ymax=mean.len+sd.len), width=.1) +
  geom_line() +
  geom_point(data=ToothGrowth, aes(x=dose, y=len, fill=supp), alpha=0.3, size=1.5) +

# Here I add the titles I want
  labs(x="Supplement Dose",
       y="Tooth Length",
       title= "The Effect of Vitamin C on Tooth Growth in Guinea Pigs",
       subtitle= "Comparing delivery methods: Orange Juice vs Ascorbic Acid",
       caption = "Based on the data from R")

Define the limits of axes

ggplot(tg.stats, aes(x=dose, y=mean.len, colour=supp)) +
  geom_point(aes(x=dose, y= mean.len), size=4, shape=18) +
  geom_errorbar(aes(ymin=mean.len-sd.len, ymax=mean.len+sd.len), width=.1) +
  geom_line() +
  geom_point(data=ToothGrowth, aes(x=dose, y=len, fill=supp), alpha=0.3, size=1.5) +

# Here I add the titles I want
  labs(x="Supplement Dose",
       y="Tooth Length",
       title= "The Effect of Vitamin C on Tooth Growth in Guinea Pigs",
       subtitle= "Comparing delivery methods: Orange Juice vs Ascorbic Acid",
       caption = "Based on the data from R") +

# Change x and y axes limits
  xlim(0, 3) +
  ylim(0, 60)

7.2 2. Change legend title

We saw in the previous session how to change legend parameters with guides() function. To change the legend title we can use the same function. See (Reference ggplot2 Guides)[http://ggplot2.tidyverse.org/reference/index.html#section-guides-axes-and-legends]

ggplot(tg.stats, aes(x=dose, y=mean.len, colour=supp)) +
  geom_point(aes(x=dose, y= mean.len), size=4, shape=18) +
  geom_errorbar(aes(ymin=mean.len-sd.len, ymax=mean.len+sd.len), width=.1) +
  geom_line() +
  geom_point(data=ToothGrowth, aes(x=dose, y=len, fill=supp), alpha=0.3, size=1.5) +

# Here I add the titles I want
  labs(x="Supplement Dose",
       y="Tooth Length",
       title= "The Effect of Vitamin C on Tooth Growth in Guinea Pigs",
       subtitle= "Comparing delivery methods: Orange Juice vs Ascorbic Acid",
       caption = "Based on the data from R") +

# Change the legend
  guides(colour=guide_legend("Supplement"), fill=guide_legend("Supplement"))

7.3 3. Jitter points

I want to customize the jitter of the points. To adjust positions look at position adjustment in ggplot2 Reference

pd <- position_dodge(0.1)

# my plot is getting very long, I'll assign it
tgp <- ggplot(tg.stats, aes(x=dose, y=mean.len, colour=supp)) +
  geom_point(aes(x=dose, y= mean.len), size=4, shape=18, position=pd) +
  geom_errorbar(aes(ymin=mean.len-sd.len, ymax=mean.len+sd.len), width=.1, position=pd) +
  geom_line(position=pd) +
  geom_point(data=ToothGrowth, aes(x=dose, y=len, fill=supp), alpha=0.3, size=1.5, position=pd) +


# Here I add the titles I want
  labs(x="Supplement Dose",
       y="Tooth Length",
       title= "The Effect of Vitamin C on Tooth Growth in Guinea Pigs",
       subtitle= "Comparing delivery methods: Orange Juice vs Ascorbic Acid",
       caption = "Based on the data from R") +

# Change the legend
  guides(colour=guide_legend("Supplement"), fill=guide_legend("Supplement"))

tgp

8.1 Customize Time Format

Now we can decide which date format we want to use for plotting. The table below explains how to provide the date format. Click here for more details.

Here some examples:

p+scale_x_date(date_labels = "%b")

p+scale_x_date(date_labels = "%Y %b %d")

p+scale_x_date(date_labels = "%m-%Y")

8.2 Customize Time Scale

To see the x-axis, let’s angle out x-axis tick labels so we can see them clearly.

p<- p + theme(axis.text.x=element_text(angle=60, hjust=1, size=8))

p + scale_x_date(date_breaks = "1 month", date_labels = "%b")

p + scale_x_date(date_breaks = "3 month", date_labels = "%b %y")

p + scale_x_date(date_breaks = "3 month", date_labels = "%b %y",
                 date_minor_breaks = "2 week") 

8.3 Limit Time Scale

You can also plot a specific period by giving a range with first and last date of the period.

p + scale_x_date(limit=c(as.Date("1994-01-01"),as.Date("1996-01-31")))